% testISS.m
%
% Given data inputs, parameter inputs, and grid will compute relevant
% growth rates (economic and ISTS) by solving a fixed point problem on the
% chain-weighted GDP growth rate.  
%
% Will then create a .txt file that can be called from latex to create a
% table.
%
% Created 10/29/2014 by Todd Messer
clear
clc
%% User Inputs
% i_y as param?
flag.iyparam = 1;
plotstuff = 0;

% Calibrated Parameters
StartYear       = [1948        ; 1988        ;1993        ; 1994];
EndYear         = [2007.50     ; 2007.50     ;2007.50     ; 2007.50];
EconGrowthRate  = [0.49965756  ;0.445451273 ;0.472815088 ;0.4635949];
InvGrowthRate   = [0.735345826 ;0.829270248 ;1.077460351 ;1.015004695];
alf             = [0.393741469 ;0.401838633 ;0.403454516 ;0.40424539];
i_y             = [0.254920809 ;0.259459451 ;0.263897188 ;0.26578];
g_y             = [0.200919601 ;0.166297844 ;0.156715936 ;0.15435684];
gammaG          = [0.057744146 ;0.124978655 ;0.209172848 ;0.197897096];

% Preliminary parameter inputs
% Parameter Values taken from model run
betin   = 0.09;             % Beta (Inverse)                            --> Estimated, Used to Calculate sigma_c

if flag.iyparam==0
    delta0  = 0.0185;
    r_ss    = 0.00625;
end

% Transform Parameters (Don't Touch)
% Use the appropriate values
bet     = (betin/100+1)^(-1);
gammaG  = gammaG/100;

% Grid to estimate gammaA
% These are the values that will be used for the experimental runs.
if flag.iyparam==1
    %     i_ygrid     = i_y;
    par     = [EconGrowthRate InvGrowthRate alf g_y gammaG i_y StartYear EndYear];

    %     par_rss    = linspace(0.005,0.009,5)';
    %     par_rss = log(1.005895932510335);
    %     par_rss = 0.00595;
    %         par_rss = 0.01;
    %     par_rss    = 0.00625;  
    par_rss = 1.015;
    
else
    alf_grid    = alf;
    d0_grid     = delta0;
    par_rss    = r_ss;
    [par_d0,par_alf] = meshgrid(d0_grid,alf_grid);
    par = [par_alf par_d0];
end

% Locations
% location to save latex file.
LatexLoc = 'C:\Users\g1txm02\Desktop\MacroFun\DSGE Model\TRY-JRC-ADD-NOMINAL-SHARES\FOMC Model 3.0\Estimation';

%% Solving for gammaA
% Loop over the parameter space
for j = 1:numel(par_rss)
    TableOut{j} = nan(0,50);
    for i = 1:size(par,1)

        % Use the current draw of parameters
        if flag.iyparam==1
            EconGrowthRate  = par(i,1);
            InvGrowthRate   = par(i,2);
            alf             = par(i,3);
            g_y             = par(i,4);
            gammaG          = par(i,5);
            i_y             = par(i,6);
            StartYear       = par(i,7);
            EndYear         = par(i,8);
            r_ss    = par_rss(j);
        else
            alf     = par(i,1);
            r_ss    = par_rss(j);
            delta0  = par(i,2);
        end

        % Solve for gammaA and gammaAx 
        if flag.iyparam==1
            gammaAdiff = @(z)fgamA_iy(z,EconGrowthRate,InvGrowthRate,g_y,i_y,gammaG);
        else
            gammaAdiff = @(z)fgamA(z,EconGrowthRate,InvGrowthRate,g_y,delta0,r_ss,alf,gammaG);
        end
        gammaA  = fzero(gammaAdiff,0.002);
        gammaAx = InvGrowthRate/100 - gammaA;

        % Calculate the other relevant variables
        if flag.iyparam==1
            delta0     = 1 - (exp(gammaAx)*((r_ss*i_y - exp(gammaA)*alf)/(i_y - alf)));
        end
        
        rk_ss     = r_ss*exp(gammaAx)-(1-delta0);
        i_k       = (1-(1-delta0)*(exp(-gammaA - gammaAx)));
        delta1    = r_ss*exp(gammaAx) - (1-delta0);
        x_y{j}(i) = alf*((1 - (1-delta0)*(exp(-gammaA-gammaAx)))/(delta1*exp(-gammaA-gammaAx)));
        i_y2      = x_y{j}(i);

        % Checks for i_y and delta0    
        if flag.iyparam==1
            if alf==i_y
                delta02{j}(i) = nan;
                zss{j}(i)  = nan;
                continue
            end
        else
            if alf==i_y2
                delta02{j}(i) = nan;
                zss{j}(i)  = nan;
                continue
            end
        end
        
        if flag.iyparam==1
%             disp(abs(i_y-i_y2))
            assert(abs(i_y-i_y2)<1e-9)
        else
            i_y = i_y2;
        end

        if flag.iyparam==0
            deltacheck  = 1-(((((i_y*exp(-gammaA-gammaAx))/(alf))*exp(r_ss + gammaAx)) - 1)/(((i_y*exp(-gammaA-gammaAx))/(alf)) - (exp(-gammaA - gammaAx))));
            assert(abs(deltacheck-delta0)<1e-9)
        end

        delta02{j}(i) = delta0;
        zss{j}(i)  = gammaA;
        
        sigmac  = (log((rk_ss + (1-delta0))*bet) - gammaAx)/gammaA;

        % Save output
        Headers = {};

        n = 0;

        TableOut{j}(i,n+1) = StartYear;             n = n+1; Headers = [Headers 'Start'];
        TableOut{j}(i,n+1) = EndYear;               n = n+1; Headers = [Headers 'End'];
        TableOut{j}(i,n+1) = EconGrowthRate;        n = n+1; Headers = [Headers '\gamma_Q'];
        TableOut{j}(i,n+1) = InvGrowthRate;         n = n+1; Headers = [Headers '\gamma_X'];
        TableOut{j}(i,n+1) = alf;                   n = n+1; Headers = [Headers 'k_y'];
        TableOut{j}(i,n+1) = i_y;                   n = n+1; Headers = [Headers 'x_y'];
        TableOut{j}(i,n+1) = g_y;                   n = n+1; Headers = [Headers 'g_y'];
        TableOut{j}(i,n+1) = 100*gammaG;            n = n+1; Headers = [Headers '\pi^g'];
%         TableOut{j}(i,n+1) = rk_ss;                 n = n+1; Headers = [Headers 'rk_ss'];
%         TableOut{j}(i,n+1) = 1/y_k;                 n = n+1; Headers = [Headers 'k_y'];
%         TableOut{j}(i,n+1) = i_k;                   n = n+1; Headers = [Headers 'i_k'];
%         TableOut{j}(i,n+1) = i_y-alf;               n = n+1; Headers = [Headers 'iy-a'];
        TableOut{j}(i,n+1) = 100*gammaA;            n = n+1; Headers = [Headers 'z_*'];
        TableOut{j}(i,n+1) = 100*gammaAx;           n = n+1; Headers = [Headers '\omega_*'];
        TableOut{j}(i,n+1) = alf*exp(gammaA)/i_y;   n = n+1; Headers = [Headers '\overline{r_*}'];
        TableOut{j}(i,n+1) = i_y^(-1)*((exp(-gammaAx)*(i_y - alf))+exp(gammaA)*alf); n=n+1; Headers = [Headers '\underline{r_*}'];
        TableOut{j}(i,n+1) = delta0;                n = n+1; Headers = [Headers '\delta_0'];
        TableOut{j}(i,n+1) = r_ss;                  n = n+1; Headers = [Headers 'r_{ss}'];
        
%         TableOut{j}(i,n+1) = alf/i_y;               n = n+1; Headers = [Headers 'alf/iy'];
%         TableOut{j}(i,n+1) = (1-exp(-gammaA-gammaAx)); n=n+1; Headers = [Headers 'Constant for alf'];
%         TableOut{j}(i,n+1) = (exp(r_ss+gammaAx) - 1); n=n+1; Headers = [Headers 'Constant for xy'];
%         TableOut{j}(i,n+1) = exp(-gammaAx)*(1-((exp(-gammaA-gammaAx)-1)*alf)/(i_y)); n=n+1; Headers = [Headers 'Real Interest Rate'];

        
%         TableOut{j}(i,n+1) = sigmac;                n = n+1; Headers = [Headers 'sigmac'];
%         TableOut{j}(i,n+1) = EconGrowthRate;        n = n+1; Headers = [Headers 'gammaQ'];
%         TableOut{j}(i,n+1) = InvGrowthRate;         n = n+1; Headers = [Headers 'gammaX'];

    end
end


% Display Table (to check)
TableOutFull = TableOut{1};
for j = 2:length(par_rss)
    TableOutFull = [TableOutFull; TableOut{j}];
end

TableDisplay = [Headers num2cell(nan(1,size(TableOutFull,2)-length(Headers))); num2cell(TableOutFull)];
disp(TableDisplay(:,cell2mat(cellfun(@(x)(all(~isnan(x))),TableDisplay(1,:),'uniformoutput',false))))

%% Ploti_y
if plotstuff
close all
colors = hsv(length(par_rss));

% \delta_0 vs. i_y, multiple r_{ss}
minz = 20; maxz = 0;
fig = figure('paperorientation','landscape','paperpositionmode','auto','position',[50 50 1200 800]);
for j = 1:length(par_rss)
    plot(par,delta02{j},'color',colors(j,:))
    leg{j} = sprintf('r_{ss} = %4.4f',par_rss(j));
    hold on
    
    minz = nanmin([minz min(zss{j})]);
    maxz = nanmax([maxz max(zss{j})]);
end
lh = line([min(par) max(par)],[0,0],'linewidth',2,'color','k');
curax = get(gca,'ylim');
vl = line([alf alf],[-100 100],'color','k','linewidth',0.5);
ylim(curax)

zssrange = sprintf('\\fontsize{14}z_{*} range: %4.4f-%4.4f',[minz maxz]);
alfrange = sprintf('\\fontsize{14}\\alpha: %3.3f',alf);
title({'Relationship Between Calibrated Parameters' zssrange alfrange},'fontsize',18)
xlabel('Investment to Output Ratio: i_y','fontsize',14)
ylabel('Steady-State Depreciation Rate:  \delta_0','fontsize',14)
xlim([min(par) max(par)])
legend(leg)
uistack(lh,'bottom')

figname = sprintf('delta vs iy - %3.3f.pdf',alf);
print(fig,'-dpdf',figname)

% \delta_0 vs. k_y, multiple r_{ss}
minz = 20; maxz = 0;
fig = figure('paperorientation','landscape','paperpositionmode','auto','position',[50 50 1200 800]);
for j = 1:length(par_rss)
    plot(k_y{j},delta02{j},'color',colors(j,:))
    leg{j} = sprintf('r_{ss} = %4.4f',par_rss(j));
    hold on
    
    minz = nanmin([minz min(zss{j})]);
    maxz = nanmax([maxz max(zss{j})]);
end
lh = line([min(k_y{j}) max(k_y{j})],[0,0],'linewidth',2,'color','k');
curax = get(gca,'ylim');
vl = line([alf alf],[-100 100],'color','k','linewidth',0.5);
ylim(curax)

zssrange = sprintf('\\fontsize{14}z_{*} range: %4.4f-%4.4f',[minz maxz]);
alfrange = sprintf('\\fontsize{14}\\alpha: %3.3f',alf);
title({'Relationship Between Calibrated Parameters' zssrange alfrange},'fontsize',18)
xlabel('Capital to Output Ratio: k_y','fontsize',14)
ylabel('Steady-State Depreciation Rate:  \delta_0','fontsize',14)
xlim([min(k_y{j}) max(k_y{j})])
legend(leg)
uistack(lh,'bottom')

figname = sprintf('delta vs ky - %3.3f.pdf',alf);
print(fig,'-dpdf',figname)

% i_y vs. k_y, multiple r_{ss}
minz = 20; maxz = 0;
fig = figure('paperorientation','landscape','paperpositionmode','auto','position',[50 50 1200 800]);
for j = 1:length(par_rss)
    plot(k_y{j},par,'color',colors(j,:))
    leg{j} = sprintf('r_{ss} = %4.4f',par_rss(j));
    hold on
    
    minz = nanmin([minz min(zss{j})]);
    maxz = nanmax([maxz max(zss{j})]);
end
zssrange = sprintf('\\fontsize{14}z_{*} range: %4.4f-%4.4f',[minz maxz]);
alfrange = sprintf('\\fontsize{14}\\alpha: %3.3f',alf);
title({'Relationship Between Calibrated Parameters' zssrange alfrange},'fontsize',18)
xlabel('Capital to Output Ratio: k_y','fontsize',14)
ylabel('Investment to Output Ratio: i_y','fontsize',14)
xlim([min(k_y{j}) max(k_y{j})])
legend(leg)
uistack(lh,'bottom')

figname = sprintf('i_y vs ky - %3.3f.pdf',alf);
print(fig,'-dpdf',figname)

end

%% Create Latex file
curcd = cd;
cd(LatexLoc)
fid = fopen('ParameterOutputLatex.txt','w');

% Headers
fprintf(fid,'\\begin{tabular}{%s} \\toprule\r\n',repmat('l',size(Headers)));
fprintf(fid,'$%s$',Headers{1});
for i = 2:length(Headers)
    fprintf(fid,'& $%s$',Headers{i});
end
fprintf(fid,'\\\\ \\midrule \r\n');

% Values
for i = 1:size(par,1)
    fprintf(fid,'%5.0f',TableOutFull(i,1));
    fprintf(fid,'& %5.1f',TableOutFull(i,2));
    
    for j = 3:length(Headers)
        fprintf(fid,'& %5.4f',TableOutFull(i,j));
    end
    fprintf(fid,'\\\\ \r\n');
end
fprintf(fid,'\\bottomrule \\\\ \r\n \\end{tabular}');

fclose(fid);
cd(curcd);